Libraries and paths

# Libraries
library(easypackages)
libraries("here","tidyverse","ggeasy","patchwork","reshape2")

fontSize = 16
dotSize = 2
alphaVal = 0.5
dot_cols2use = c("#404040","#a3142e")
dtype_cols2use = c("#7ab5ed","#171c9e")

# Paths
codedir = here("code")
plotdir = here("plots")
datadir = here("data")


# functions that will be used over and over again

# make_scatterplot
make_scatterplot <- function(data2plot, x_var, y_var, strat_var, 
                             xlabel2use, ylabel2use, title2use,
                             vlineval=NULL,plot_fit_line=FALSE,plot_diag_line=FALSE,
                             dotSize=2, fontSize=16, alphaVal=0.5,
                             cols2use=c("#404040","#a3142e")
                             ){
  
  p = ggplot(data = data2plot, aes_string(x = x_var, y = y_var, colour= strat_var))
  
  if (is.null(alphaVal)){
    p = p + geom_point(size = dotSize)
  } else{
    p = p + geom_point(size = dotSize, alpha = alphaVal) 
  }
  
  p = p + 
    scale_colour_manual(values=cols2use) + 
    xlab(xlabel2use) +  
    ylab(ylabel2use) + 
    ggtitle(title2use) +
    easy_center_title() +
    guides(colour="none")+
    theme_bw() + 
    theme(plot.title = element_text(hjust = 0.5, size=fontSize),
        strip.text.y = element_text(size=fontSize),
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize-2),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text = element_text(size = fontSize))
  
  if (!is.null(vlineval)){
    p = p + geom_vline(xintercept = vlineval) 
  }
  
  if (plot_fit_line){
    p = p + geom_smooth(method = lm, se=FALSE, aes_string(group=strat_var))
  }
  
  if (plot_diag_line){
    p = p + geom_abline(intercept = 0, slope = 1, colour = "black")
  }
  
  return(p)
  
} # function make_scatterplot

# out_of_sample_prediction
out_of_sample_prediction <- function(trn_data, val_data, dv_var, pred_var){
  
  form2use = as.formula(sprintf("%s ~ %s", dv_var, pred_var))
  mod2use = lm(formula = form2use, data = trn_data)
  val_data$pred = predict(mod2use, val_data)
  
  # compute R-squared
  rss = sum((val_data[,"pred"] - val_data[,dv_var])^2)  ## residual sum of squares
  tss = sum((val_data[,dv_var] - mean(val_data[,dv_var]))^2)  ## total sum of squares
  rsq = 1 - (rss/tss)
  
  res_list = list(trn_data = trn_data,
                  val_data = val_data,
                  mod2use = mod2use, 
                  rsq = rsq)
  return(res_list)
} # out_of_sample_prediction

Read in data

g_eeg_H_data = read.csv(file.path(datadir,"tidy_H_g_EEG_data.csv"))
g_lfp_H_data = read.csv(file.path(datadir,"tidy_H_g_LFP_data.csv"))
g_h_data = rbind(g_eeg_H_data, g_lfp_H_data)

g_eeg_slope_data = read.csv(file.path(datadir,"tidy_slope_g_EEG_data.csv"))
g_lfp_slope_data = read.csv(file.path(datadir,"tidy_slope_g_LFP_data.csv"))
g_slope_data = rbind(g_eeg_slope_data, g_lfp_slope_data)

g_eeg_pow_data = read.csv(file.path(datadir,"tidy_pow_g_EEG_data.csv"))
g_lfp_pow_data = read.csv(file.path(datadir,"tidy_pow_g_LFP_data.csv"))
g_pow_data = rbind(g_eeg_pow_data, g_lfp_pow_data)


camk_eeg_H_data = read.csv(file.path(datadir,"tidy_H_camk_EEG_data.csv"))
camk_lfp_H_data = read.csv(file.path(datadir,"tidy_H_camk_LFP_data.csv"))
camk_h_data = rbind(camk_eeg_H_data, camk_lfp_H_data)

camk_eeg_slope_data = read.csv(file.path(datadir,"tidy_slope_camk_EEG_data.csv"))
camk_lfp_slope_data = read.csv(file.path(datadir,"tidy_slope_camk_LFP_data.csv"))
camk_slope_data = rbind(camk_eeg_slope_data, camk_lfp_slope_data)

camk_eeg_pow_data = read.csv(file.path(datadir,"tidy_pow_camk_EEG_data.csv"))
camk_lfp_pow_data = read.csv(file.path(datadir,"tidy_pow_camk_LFP_data.csv"))
camk_pow_data = rbind(camk_eeg_pow_data, camk_lfp_pow_data)


hm4di_eeg_H_data = read.csv(file.path(datadir,"tidy_H_hm4di_EEG_data.csv"))
hm4di_lfp_H_data = read.csv(file.path(datadir,"tidy_H_hm4di_LFP_data.csv"))
hm4di_h_data = rbind(hm4di_eeg_H_data, hm4di_lfp_H_data)

hm4di_eeg_slope_data = read.csv(file.path(datadir,"tidy_slope_hm4di_EEG_data.csv"))
hm4di_lfp_slope_data = read.csv(file.path(datadir,"tidy_slope_hm4di_LFP_data.csv"))
hm4di_slope_data = rbind(hm4di_eeg_slope_data, hm4di_lfp_slope_data)

hm4di_eeg_pow_data = read.csv(file.path(datadir,"tidy_pow_hm4di_EEG_data.csv"))
hm4di_lfp_pow_data = read.csv(file.path(datadir,"tidy_pow_hm4di_LFP_data.csv"))
hm4di_pow_data = rbind(hm4di_eeg_pow_data, hm4di_lfp_pow_data)

Plot Firing Rate by g

# Scatterplot of Firing rate by g
data2plot = g_h_data %>%  
  filter(dset=="training") %>% 
  filter(measure=="EEG") %>%
  select(fr_exc, fr_inh, param) %>% 
  distinct() %>% 
  melt(id.vars=c("param"))

vlineval = 0.0885
yLabel = "Firing Rate (Hz)"
xLabel = expression("g =" ~ g[E]/g[I])
title2use = xLabel
cols2use = dot_cols2use

FR_g_plot = make_scatterplot(data2plot = data2plot, 
                            x_var = "param",
                            y_var = "value",
                            strat_var = "variable", 
                            vlineval = vlineval, 
                            xlabel2use = xLabel, 
                            ylabel2use = yLabel, 
                            title2use = title2use, 
                            cols2use = cols2use, 
                            alphaVal = NULL)
FR_g_plot

# correlation with g in E neurons
data2use = data2plot %>% filter(variable=="fr_exc")
res = cor.test(data2use[,"param"], data2use[,"value"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use[, "param"] and data2use[, "value"]
## t = 19.051, df = 68, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8704428 0.9482165
## sample estimates:
##       cor 
## 0.9177152
res$p.value
## [1] 5.698427e-29
# correlation with g in I neurons
data2use = data2plot %>% filter(variable=="fr_inh")
res = cor.test(data2use[,"param"], data2use[,"value"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use[, "param"] and data2use[, "value"]
## t = 67.452, df = 68, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9880972 0.9954156
## sample estimates:
##       cor 
## 0.9926099
res$p.value
## [1] 5.032159e-64

Plot Firing Rate by Camk El

# Scatterplot of Firing rate by Camk El
data2plot = camk_h_data %>%  
  filter(dset=="training") %>% 
  filter(measure=="EEG") %>%
  select(fr_exc, fr_inh, param) %>% 
  distinct() %>% 
  melt(id.vars=c("param"))

vlineval = -70
yLabel = "Firing Rate (Hz)"
xLabel = "Resting Potential (mV)"
title2use = "CamkII-hM3D(Gq)"
cols2use = dot_cols2use

FR_camk_El_plot = make_scatterplot(data2plot = data2plot, 
                            x_var = "param",
                            y_var = "value",
                            strat_var = "variable", 
                            vlineval = vlineval, 
                            xlabel2use = xLabel, 
                            ylabel2use = yLabel, 
                            title2use = title2use, 
                            cols2use = cols2use,
                            alphaVal = NULL)
FR_camk_El_plot

# correlation with g in E neurons
data2use = data2plot %>% filter(variable=="fr_exc")
res = cor.test(data2use[,"param"], data2use[,"value"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use[, "param"] and data2use[, "value"]
## t = 52.206, df = 49, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.984434 0.994953
## sample estimates:
##       cor 
## 0.9911301
res$p.value
## [1] 1.289149e-44
# correlation with g in I neurons
data2use = data2plot %>% filter(variable=="fr_inh")
res = cor.test(data2use[,"param"], data2use[,"value"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use[, "param"] and data2use[, "value"]
## t = 59.298, df = 49, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9878897 0.9960781
## sample estimates:
##       cor 
## 0.9931044
res$p.value
## [1] 2.760435e-47

Plot Firing Rate by hM4Di El

# Scatterplot of Firing rate by hM4Di El
data2plot = hm4di_h_data %>%  
  filter(dset=="training") %>% 
  filter(measure=="EEG") %>%
  select(fr_exc, fr_inh, param) %>% 
  distinct() %>% 
  melt(id.vars=c("param"))

vlineval = -70
yLabel = "Firing Rate (Hz)"
xLabel = "Resting Potential (mV)"
title2use = "hM4Di"
cols2use = dot_cols2use

FR_hm4di_El_plot = make_scatterplot(data2plot = data2plot, 
                            x_var = "param",
                            y_var = "value",
                            strat_var = "variable", 
                            vlineval = vlineval, 
                            xlabel2use = xLabel, 
                            ylabel2use = yLabel, 
                            title2use = title2use, 
                            cols2use = cols2use,
                            alphaVal = NULL)
FR_hm4di_El_plot

# correlation with g in E neurons
data2use = data2plot %>% filter(variable=="fr_exc")
res = cor.test(data2use[,"param"], data2use[,"value"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use[, "param"] and data2use[, "value"]
## t = 5.5862, df = 9, p-value = 0.0003402
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5962093 0.9688472
## sample estimates:
##       cor 
## 0.8809954
res$p.value
## [1] 0.0003402203
# correlation with g in I neurons
data2use = data2plot %>% filter(variable=="fr_inh")
res = cor.test(data2use[,"param"], data2use[,"value"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use[, "param"] and data2use[, "value"]
## t = 17.813, df = 9, p-value = 2.515e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9456054 0.9965087
## sample estimates:
##       cor 
## 0.9861128
res$p.value
## [1] 2.5149e-08

Analysis of E:I ratio simulation

# ==============================================================================
# H
# Scatterplot against g
data2plot = g_h_data %>% filter(dset=="training")

vlineval = 0.0885
yLabel = "H"
xLabel = expression("g =" ~ g[E]/g[I])
title2use = xLabel
cols2use = dtype_cols2use
plot_fit_line = TRUE

H_g_plot = make_scatterplot(data2plot = data2plot, 
                            x_var = "param",
                            y_var = "value",
                            strat_var = "measure", 
                            vlineval = vlineval, 
                            xlabel2use = xLabel, 
                            ylabel2use = yLabel, 
                            title2use = title2use, 
                            cols2use = cols2use, 
                            plot_fit_line = plot_fit_line)
H_g_plot

# correlation with g in LFP
data2use = data2plot %>% filter(measure=="LFP")
res = cor.test(data2use[,"param"], data2use[,"value"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use[, "param"] and data2use[, "value"]
## t = -11.035, df = 68, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.8719076 -0.6973305
## sample estimates:
##       cor 
## -0.801049
res$p.value
## [1] 8.36119e-17
# correlation with g in EEG
data2use = data2plot %>% filter(measure=="EEG")
res = cor.test(data2use[,"param"], data2use[,"value"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use[, "param"] and data2use[, "value"]
## t = -11.365, df = 68, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.8774729 -0.7093028
## sample estimates:
##        cor 
## -0.8093794
res$p.value
## [1] 2.262083e-17
# ==============================================================================
# out of sample prediction
trn_eeg_data = g_h_data %>% filter(dset=="training") %>% filter(measure=="EEG")
val_eeg_data = g_h_data %>% filter(dset=="validation") %>% filter(measure=="EEG")
oos_res_eeg = out_of_sample_prediction(trn_data = trn_eeg_data, 
                                   val_data = val_eeg_data, 
                                   dv_var="fr_tot", 
                                   pred_var="value")
print(sprintf("EEG R-sq = %f",oos_res_eeg$rsq))
## [1] "EEG R-sq = 0.558346"
trn_lfp_data = g_h_data %>% filter(dset=="training") %>% filter(measure=="LFP")
val_lfp_data = g_h_data %>% filter(dset=="validation") %>% filter(measure=="LFP")
oos_res_lfp = out_of_sample_prediction(trn_data = trn_lfp_data, 
                                   val_data = val_lfp_data, 
                                   dv_var="fr_tot", 
                                   pred_var="value")
print(sprintf("LFP R-sq = %f",oos_res_lfp$rsq))
## [1] "LFP R-sq = 0.671406"
# make plot of actual vs predicted
data2plot = rbind(oos_res_eeg$val_data, oos_res_lfp$val_data)

vlineval = NULL
yLabel = "Predicted Firing Rate (Hz)"
xLabel = "Actual Firing Rate (Hz)"
title2use = expression("g =" ~ g[E]/g[I])
cols2use = dtype_cols2use
plot_diag_line = TRUE

H_g_oos_plot = make_scatterplot(data2plot = data2plot, 
                            x_var = "fr_tot",
                            y_var = "pred",
                            strat_var = "measure", 
                            vlineval = vlineval, 
                            xlabel2use = xLabel, 
                            ylabel2use = yLabel, 
                            title2use = title2use, 
                            cols2use = cols2use, 
                            plot_diag_line = plot_diag_line)
H_g_oos_plot

# ==============================================================================
# 1/f slope
# Scatterplot against g
data2plot = g_slope_data %>% filter(dset=="training")

vlineval = 0.0885
yLabel = "1/f slope"
xLabel = expression("g =" ~ g[E]/g[I])
title2use = xLabel
cols2use = dtype_cols2use
plot_fit_line = TRUE

slope_g_plot = make_scatterplot(data2plot = data2plot, 
                            x_var = "param",
                            y_var = "value",
                            strat_var = "measure", 
                            vlineval = vlineval, 
                            xlabel2use = xLabel, 
                            ylabel2use = yLabel, 
                            title2use = title2use, 
                            cols2use = cols2use, 
                            plot_fit_line = plot_fit_line)
slope_g_plot

# correlation with g in LFP
data2use = data2plot %>% filter(measure=="LFP")
res = cor.test(data2use[,"param"], data2use[,"value"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use[, "param"] and data2use[, "value"]
## t = 16.573, df = 68, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8362499 0.9338236
## sample estimates:
##       cor 
## 0.8953001
res$p.value
## [1] 1.411939e-25
# correlation with g in EEG
data2use = data2plot %>% filter(measure=="EEG")
res = cor.test(data2use[,"param"], data2use[,"value"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use[, "param"] and data2use[, "value"]
## t = 17.487, df = 68, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8501953 0.9397321
## sample estimates:
##       cor 
## 0.9044783
res$p.value
## [1] 7.280804e-27
# ==============================================================================
# out of sample prediction
trn_eeg_data = g_slope_data %>% filter(dset=="training") %>% filter(measure=="EEG")
val_eeg_data = g_slope_data %>% filter(dset=="validation") %>% filter(measure=="EEG")
oos_res_eeg = out_of_sample_prediction(trn_data = trn_eeg_data, 
                                   val_data = val_eeg_data, 
                                   dv_var="fr_tot", 
                                   pred_var="value")
print(sprintf("EEG R-sq = %f",oos_res_eeg$rsq))
## [1] "EEG R-sq = 0.747769"
trn_lfp_data = g_slope_data %>% filter(dset=="training") %>% filter(measure=="LFP")
val_lfp_data = g_slope_data %>% filter(dset=="validation") %>% filter(measure=="LFP")
oos_res_lfp = out_of_sample_prediction(trn_data = trn_lfp_data, 
                                   val_data = val_lfp_data, 
                                   dv_var="fr_tot", 
                                   pred_var="value")
print(sprintf("LFP R-sq = %f",oos_res_lfp$rsq))
## [1] "LFP R-sq = 0.711435"
# make plot of actual vs predicted
data2plot = rbind(oos_res_eeg$val_data, oos_res_lfp$val_data)

vlineval = NULL
yLabel = "Predicted Firing Rate (Hz)"
xLabel = "Actual Firing Rate (Hz)"
title2use = expression("g =" ~ g[E]/g[I])
cols2use = dtype_cols2use
plot_diag_line = TRUE

slope_g_oos_plot = make_scatterplot(data2plot = data2plot, 
                            x_var = "fr_tot",
                            y_var = "pred",
                            strat_var = "measure", 
                            vlineval = vlineval, 
                            xlabel2use = xLabel, 
                            ylabel2use = yLabel, 
                            title2use = title2use, 
                            cols2use = cols2use, 
                            plot_diag_line = plot_diag_line)
slope_g_oos_plot

# ==============================================================================
# Power
# Scatterplot against g
data2plot = g_pow_data %>% filter(dset=="training")

vlineval = 0.0885
yLabel = "Broadband Power (a.u.)"
xLabel = expression("g =" ~ g[E]/g[I])
title2use = xLabel
cols2use = dtype_cols2use
plot_fit_line = TRUE

pow_g_plot = make_scatterplot(data2plot = data2plot, 
                            x_var = "param",
                            y_var = "value",
                            strat_var = "measure", 
                            vlineval = vlineval, 
                            xlabel2use = xLabel, 
                            ylabel2use = yLabel, 
                            title2use = title2use, 
                            cols2use = cols2use, 
                            plot_fit_line = plot_fit_line)
pow_g_plot

# correlation with g in LFP
data2use = data2plot %>% filter(measure=="LFP")
res = cor.test(data2use[,"param"], data2use[,"value"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use[, "param"] and data2use[, "value"]
## t = 34.722, df = 68, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9566735 0.9831489
## sample estimates:
##       cor 
## 0.9729376
res$p.value
## [1] 5.370013e-45
# correlation with g in EEG
data2use = data2plot %>% filter(measure=="EEG")
res = cor.test(data2use[,"param"], data2use[,"value"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use[, "param"] and data2use[, "value"]
## t = 34.216, df = 68, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9554483 0.9826658
## sample estimates:
##       cor 
## 0.9721658
res$p.value
## [1] 1.37964e-44
# ==============================================================================
# out of sample prediction
trn_eeg_data = g_pow_data %>% filter(dset=="training") %>% filter(measure=="EEG")
val_eeg_data = g_pow_data %>% filter(dset=="validation") %>% filter(measure=="EEG")
oos_res_eeg = out_of_sample_prediction(trn_data = trn_eeg_data, 
                                   val_data = val_eeg_data, 
                                   dv_var="fr_tot", 
                                   pred_var="value")
print(sprintf("EEG R-sq = %f",oos_res_eeg$rsq))
## [1] "EEG R-sq = 0.911498"
trn_lfp_data = g_pow_data %>% filter(dset=="training") %>% filter(measure=="LFP")
val_lfp_data = g_pow_data %>% filter(dset=="validation") %>% filter(measure=="LFP")
oos_res_lfp = out_of_sample_prediction(trn_data = trn_lfp_data, 
                                   val_data = val_lfp_data, 
                                   dv_var="fr_tot", 
                                   pred_var="value")
print(sprintf("LFP R-sq = %f",oos_res_lfp$rsq))
## [1] "LFP R-sq = 0.912992"
# make plot of actual vs predicted
data2plot = rbind(oos_res_eeg$val_data, oos_res_lfp$val_data)

vlineval = NULL
yLabel = "Predicted Firing Rate (Hz)"
xLabel = "Actual Firing Rate (Hz)"
title2use = expression("g =" ~ g[E]/g[I])
cols2use = dtype_cols2use
plot_diag_line = TRUE

pow_g_oos_plot = make_scatterplot(data2plot = data2plot, 
                            x_var = "fr_tot",
                            y_var = "pred",
                            strat_var = "measure", 
                            vlineval = vlineval, 
                            xlabel2use = xLabel, 
                            ylabel2use = yLabel, 
                            title2use = title2use, 
                            cols2use = cols2use, 
                            plot_diag_line = plot_diag_line)
pow_g_oos_plot

Analysis of CamK

title2use = "CamkII-hM3D(Gq)"

# ==============================================================================
# H
# Scatterplot against El
data2plot = camk_h_data %>% filter(dset=="training")

vlineval = -70
yLabel = "H"
xLabel = "Resting Potential (mV)"
cols2use = dtype_cols2use
plot_fit_line = TRUE

H_camk_El_plot = make_scatterplot(data2plot = data2plot, 
                            x_var = "param",
                            y_var = "value",
                            strat_var = "measure", 
                            vlineval = vlineval, 
                            xlabel2use = xLabel, 
                            ylabel2use = yLabel, 
                            title2use = title2use, 
                            cols2use = cols2use, 
                            plot_fit_line = plot_fit_line)
H_camk_El_plot

# correlation with El in LFP
data2use = data2plot %>% filter(measure=="LFP")
res = cor.test(data2use[,"param"], data2use[,"value"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use[, "param"] and data2use[, "value"]
## t = -37.543, df = 49, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.9903432 -0.9703590
## sample estimates:
##        cor 
## -0.9830581
res$p.value
## [1] 9.037516e-38
# correlation with El in EEG
data2use = data2plot %>% filter(measure=="EEG")
res = cor.test(data2use[,"param"], data2use[,"value"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use[, "param"] and data2use[, "value"]
## t = -47.763, df = 49, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.9939838 -0.9814634
## sample estimates:
##        cor 
## -0.9894305
res$p.value
## [1] 9.272795e-43
# ==============================================================================
# out of sample prediction
trn_eeg_data = camk_h_data %>% filter(dset=="training") %>% filter(measure=="EEG")
val_eeg_data = camk_h_data %>% filter(dset=="validation") %>% filter(measure=="EEG")
oos_res_eeg = out_of_sample_prediction(trn_data = trn_eeg_data, 
                                   val_data = val_eeg_data, 
                                   dv_var="fr_tot", 
                                   pred_var="value")
print(sprintf("EEG R-sq = %f",oos_res_eeg$rsq))
## [1] "EEG R-sq = 0.953605"
trn_lfp_data = camk_h_data %>% filter(dset=="training") %>% filter(measure=="LFP")
val_lfp_data = camk_h_data %>% filter(dset=="validation") %>% filter(measure=="LFP")
oos_res_lfp = out_of_sample_prediction(trn_data = trn_lfp_data, 
                                   val_data = val_lfp_data, 
                                   dv_var="fr_tot", 
                                   pred_var="value")
print(sprintf("LFP R-sq = %f",oos_res_lfp$rsq))
## [1] "LFP R-sq = 0.952495"
# make plot of actual vs predicted
data2plot = rbind(oos_res_eeg$val_data, oos_res_lfp$val_data)

vlineval = NULL
yLabel = "Predicted Firing Rate (Hz)"
xLabel = "Actual Firing Rate (Hz)"
cols2use = dtype_cols2use
plot_diag_line = TRUE

H_camk_El_oos_plot = make_scatterplot(data2plot = data2plot, 
                            x_var = "fr_tot",
                            y_var = "pred",
                            strat_var = "measure", 
                            vlineval = vlineval, 
                            xlabel2use = xLabel, 
                            ylabel2use = yLabel, 
                            title2use = title2use, 
                            cols2use = cols2use, 
                            plot_diag_line = plot_diag_line)
H_camk_El_oos_plot

# ==============================================================================
# 1/f slope
# Scatterplot against El
data2plot = camk_slope_data %>% filter(dset=="training")

vlineval = -70
yLabel = "1/f slope"
xLabel = "Resting Potential (mV)"
cols2use = dtype_cols2use
plot_fit_line = TRUE

slope_camk_El_plot = make_scatterplot(data2plot = data2plot, 
                            x_var = "param",
                            y_var = "value",
                            strat_var = "measure", 
                            vlineval = vlineval, 
                            xlabel2use = xLabel, 
                            ylabel2use = yLabel, 
                            title2use = title2use, 
                            cols2use = cols2use, 
                            plot_fit_line = plot_fit_line)
slope_camk_El_plot

# correlation with g in LFP
data2use = data2plot %>% filter(measure=="LFP")
res = cor.test(data2use[,"param"], data2use[,"value"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use[, "param"] and data2use[, "value"]
## t = 12.396, df = 49, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7830842 0.9244918
## sample estimates:
##       cor 
## 0.8707545
res$p.value
## [1] 1.014247e-16
# correlation with g in EEG
data2use = data2plot %>% filter(measure=="EEG")
res = cor.test(data2use[,"param"], data2use[,"value"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use[, "param"] and data2use[, "value"]
## t = 13.587, df = 49, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8123856 0.9353838
## sample estimates:
##       cor 
## 0.8889508
res$p.value
## [1] 3.062337e-18
# ==============================================================================
# out of sample prediction
trn_eeg_data = camk_slope_data %>% filter(dset=="training") %>% filter(measure=="EEG")
val_eeg_data = camk_slope_data %>% filter(dset=="validation") %>% filter(measure=="EEG")
oos_res_eeg = out_of_sample_prediction(trn_data = trn_eeg_data, 
                                   val_data = val_eeg_data, 
                                   dv_var="fr_tot", 
                                   pred_var="value")
print(sprintf("EEG R-sq = %f",oos_res_eeg$rsq))
## [1] "EEG R-sq = 0.864963"
trn_lfp_data = camk_slope_data %>% filter(dset=="training") %>% filter(measure=="LFP")
val_lfp_data = camk_slope_data %>% filter(dset=="validation") %>% filter(measure=="LFP")
oos_res_lfp = out_of_sample_prediction(trn_data = trn_lfp_data, 
                                   val_data = val_lfp_data, 
                                   dv_var="fr_tot", 
                                   pred_var="value")
print(sprintf("LFP R-sq = %f",oos_res_lfp$rsq))
## [1] "LFP R-sq = 0.849740"
# make plot of actual vs predicted
data2plot = rbind(oos_res_eeg$val_data, oos_res_lfp$val_data)

vlineval = NULL
yLabel = "Predicted Firing Rate (Hz)"
xLabel = "Actual Firing Rate (Hz)"
cols2use = dtype_cols2use
plot_diag_line = TRUE

slope_camk_El_oos_plot = make_scatterplot(data2plot = data2plot, 
                            x_var = "fr_tot",
                            y_var = "pred",
                            strat_var = "measure", 
                            vlineval = vlineval, 
                            xlabel2use = xLabel, 
                            ylabel2use = yLabel, 
                            title2use = title2use, 
                            cols2use = cols2use, 
                            plot_diag_line = plot_diag_line)
slope_camk_El_oos_plot

# ==============================================================================
# Power
# Scatterplot against g
data2plot = camk_pow_data %>% filter(dset=="training")

vlineval = -70
yLabel = "Broadband Power (a.u.)"
xLabel = "Resting Potential (mV)"
cols2use = dtype_cols2use
plot_fit_line = TRUE

pow_camk_El_plot = make_scatterplot(data2plot = data2plot, 
                            x_var = "param",
                            y_var = "value",
                            strat_var = "measure", 
                            vlineval = vlineval, 
                            xlabel2use = xLabel, 
                            ylabel2use = yLabel, 
                            title2use = title2use, 
                            cols2use = cols2use, 
                            plot_fit_line = plot_fit_line)
pow_camk_El_plot

# correlation with g in LFP
data2use = data2plot %>% filter(measure=="LFP")
res = cor.test(data2use[,"param"], data2use[,"value"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use[, "param"] and data2use[, "value"]
## t = 50.2, df = 49, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9831882 0.9945468
## sample estimates:
##       cor 
## 0.9904176
res$p.value
## [1] 8.490014e-44
# correlation with g in EEG
data2use = data2plot %>% filter(measure=="EEG")
res = cor.test(data2use[,"param"], data2use[,"value"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use[, "param"] and data2use[, "value"]
## t = 48.374, df = 49, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9819199 0.9941328
## sample estimates:
##       cor 
## 0.9896919
res$p.value
## [1] 5.036442e-43
# ==============================================================================
# out of sample prediction
trn_eeg_data = camk_pow_data %>% filter(dset=="training") %>% filter(measure=="EEG")
val_eeg_data = camk_pow_data %>% filter(dset=="validation") %>% filter(measure=="EEG")
oos_res_eeg = out_of_sample_prediction(trn_data = trn_eeg_data, 
                                   val_data = val_eeg_data, 
                                   dv_var="fr_tot", 
                                   pred_var="value")
print(sprintf("EEG R-sq = %f",oos_res_eeg$rsq))
## [1] "EEG R-sq = 0.969967"
trn_lfp_data = camk_pow_data %>% filter(dset=="training") %>% filter(measure=="LFP")
val_lfp_data = camk_pow_data %>% filter(dset=="validation") %>% filter(measure=="LFP")
oos_res_lfp = out_of_sample_prediction(trn_data = trn_lfp_data, 
                                   val_data = val_lfp_data, 
                                   dv_var="fr_tot", 
                                   pred_var="value")
print(sprintf("LFP R-sq = %f",oos_res_lfp$rsq))
## [1] "LFP R-sq = 0.972305"
# make plot of actual vs predicted
data2plot = rbind(oos_res_eeg$val_data, oos_res_lfp$val_data)

vlineval = NULL
yLabel = "Predicted Firing Rate (Hz)"
xLabel = "Actual Firing Rate (Hz)"
cols2use = dtype_cols2use
plot_diag_line = TRUE

pow_camk_El_oos_plot = make_scatterplot(data2plot = data2plot, 
                            x_var = "fr_tot",
                            y_var = "pred",
                            strat_var = "measure", 
                            vlineval = vlineval, 
                            xlabel2use = xLabel, 
                            ylabel2use = yLabel, 
                            title2use = title2use, 
                            cols2use = cols2use, 
                            plot_diag_line = plot_diag_line)
pow_camk_El_oos_plot

Analysis of hM4Di

title2use = "hM4Di"

# ==============================================================================
# H
# Scatterplot against El
data2plot = hm4di_h_data %>% filter(dset=="training")

vlineval = -70
yLabel = "H"
xLabel = "Resting Potential (mV)"
cols2use = dtype_cols2use
plot_fit_line = TRUE

H_hm4di_El_plot = make_scatterplot(data2plot = data2plot, 
                            x_var = "param",
                            y_var = "value",
                            strat_var = "measure", 
                            vlineval = vlineval, 
                            xlabel2use = xLabel, 
                            ylabel2use = yLabel, 
                            title2use = title2use, 
                            cols2use = cols2use, 
                            plot_fit_line = plot_fit_line)
H_hm4di_El_plot

# correlation with El in LFP
data2use = data2plot %>% filter(measure=="LFP")
res = cor.test(data2use[,"param"], data2use[,"value"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use[, "param"] and data2use[, "value"]
## t = -10.598, df = 9, p-value = 2.203e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.9904076 -0.8569264
## sample estimates:
##        cor 
## -0.9621892
res$p.value
## [1] 2.202948e-06
# correlation with El in EEG
data2use = data2plot %>% filter(measure=="EEG")
res = cor.test(data2use[,"param"], data2use[,"value"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use[, "param"] and data2use[, "value"]
## t = -25.916, df = 9, p-value = 9.144e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.9983368 -0.9737372
## sample estimates:
##        cor 
## -0.9933664
res$p.value
## [1] 9.144413e-10
# ==============================================================================
# out of sample prediction
trn_eeg_data = hm4di_h_data %>% filter(dset=="training") %>% filter(measure=="EEG")
val_eeg_data = hm4di_h_data %>% filter(dset=="validation") %>% filter(measure=="EEG")
oos_res_eeg = out_of_sample_prediction(trn_data = trn_eeg_data, 
                                   val_data = val_eeg_data, 
                                   dv_var="fr_tot", 
                                   pred_var="value")
print(sprintf("EEG R-sq = %f",oos_res_eeg$rsq))
## [1] "EEG R-sq = 0.919404"
trn_lfp_data = hm4di_h_data %>% filter(dset=="training") %>% filter(measure=="LFP")
val_lfp_data = hm4di_h_data %>% filter(dset=="validation") %>% filter(measure=="LFP")
oos_res_lfp = out_of_sample_prediction(trn_data = trn_lfp_data, 
                                   val_data = val_lfp_data, 
                                   dv_var="fr_tot", 
                                   pred_var="value")
print(sprintf("LFP R-sq = %f",oos_res_lfp$rsq))
## [1] "LFP R-sq = 0.947733"
# make plot of actual vs predicted
data2plot = rbind(oos_res_eeg$val_data, oos_res_lfp$val_data)

vlineval = NULL
yLabel = "Predicted Firing Rate (Hz)"
xLabel = "Actual Firing Rate (Hz)"
cols2use = dtype_cols2use
plot_diag_line = TRUE

H_hm4di_El_oos_plot = make_scatterplot(data2plot = data2plot, 
                            x_var = "fr_tot",
                            y_var = "pred",
                            strat_var = "measure", 
                            vlineval = vlineval, 
                            xlabel2use = xLabel, 
                            ylabel2use = yLabel, 
                            title2use = title2use, 
                            cols2use = cols2use, 
                            plot_diag_line = plot_diag_line)
H_hm4di_El_oos_plot

# ==============================================================================
# 1/f slope
# Scatterplot against El
data2plot = hm4di_slope_data %>% filter(dset=="training")

vlineval = -70
yLabel = "1/f slope"
xLabel = "Resting Potential (mV)"
cols2use = dtype_cols2use
plot_fit_line = TRUE

slope_hm4di_El_plot = make_scatterplot(data2plot = data2plot, 
                            x_var = "param",
                            y_var = "value",
                            strat_var = "measure", 
                            vlineval = vlineval, 
                            xlabel2use = xLabel, 
                            ylabel2use = yLabel, 
                            title2use = title2use, 
                            cols2use = cols2use, 
                            plot_fit_line = plot_fit_line)
slope_hm4di_El_plot

# correlation with g in LFP
data2use = data2plot %>% filter(measure=="LFP")
res = cor.test(data2use[,"param"], data2use[,"value"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use[, "param"] and data2use[, "value"]
## t = 4.5528, df = 9, p-value = 0.00138
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4711351 0.9560171
## sample estimates:
##       cor 
## 0.8350161
res$p.value
## [1] 0.001380469
# correlation with g in EEG
data2use = data2plot %>% filter(measure=="EEG")
res = cor.test(data2use[,"param"], data2use[,"value"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use[, "param"] and data2use[, "value"]
## t = 4.3696, df = 9, p-value = 0.001798
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4442079 0.9529888
## sample estimates:
##       cor 
## 0.8244045
res$p.value
## [1] 0.001798067
# ==============================================================================
# out of sample prediction
trn_eeg_data = hm4di_slope_data %>% filter(dset=="training") %>% filter(measure=="EEG")
val_eeg_data = hm4di_slope_data %>% filter(dset=="validation") %>% filter(measure=="EEG")
oos_res_eeg = out_of_sample_prediction(trn_data = trn_eeg_data, 
                                   val_data = val_eeg_data, 
                                   dv_var="fr_tot", 
                                   pred_var="value")
print(sprintf("EEG R-sq = %f",oos_res_eeg$rsq))
## [1] "EEG R-sq = 0.658198"
trn_lfp_data = hm4di_slope_data %>% filter(dset=="training") %>% filter(measure=="LFP")
val_lfp_data = hm4di_slope_data %>% filter(dset=="validation") %>% filter(measure=="LFP")
oos_res_lfp = out_of_sample_prediction(trn_data = trn_lfp_data, 
                                   val_data = val_lfp_data, 
                                   dv_var="fr_tot", 
                                   pred_var="value")
print(sprintf("LFP R-sq = %f",oos_res_lfp$rsq))
## [1] "LFP R-sq = 0.676430"
# make plot of actual vs predicted
data2plot = rbind(oos_res_eeg$val_data, oos_res_lfp$val_data)

vlineval = NULL
yLabel = "Predicted Firing Rate (Hz)"
xLabel = "Actual Firing Rate (Hz)"
cols2use = dtype_cols2use
plot_diag_line = TRUE

slope_hm4di_El_oos_plot = make_scatterplot(data2plot = data2plot, 
                            x_var = "fr_tot",
                            y_var = "pred",
                            strat_var = "measure", 
                            vlineval = vlineval, 
                            xlabel2use = xLabel, 
                            ylabel2use = yLabel, 
                            title2use = title2use, 
                            cols2use = cols2use, 
                            plot_diag_line = plot_diag_line)
slope_hm4di_El_oos_plot

# ==============================================================================
# Power
# Scatterplot against g
data2plot = hm4di_pow_data %>% filter(dset=="training")

vlineval = -70
yLabel = "Broadband Power (a.u.)"
xLabel = "Resting Potential (mV)"
cols2use = dtype_cols2use
plot_fit_line = TRUE

pow_hm4di_El_plot = make_scatterplot(data2plot = data2plot, 
                            x_var = "param",
                            y_var = "value",
                            strat_var = "measure", 
                            vlineval = vlineval, 
                            xlabel2use = xLabel, 
                            ylabel2use = yLabel, 
                            title2use = title2use, 
                            cols2use = cols2use, 
                            plot_fit_line = plot_fit_line)
pow_hm4di_El_plot

# correlation with g in LFP
data2use = data2plot %>% filter(measure=="LFP")
res = cor.test(data2use[,"param"], data2use[,"value"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use[, "param"] and data2use[, "value"]
## t = 6.5724, df = 9, p-value = 0.0001025
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6820298 0.9766279
## sample estimates:
##       cor 
## 0.9097119
res$p.value
## [1] 0.0001024801
# correlation with g in EEG
data2use = data2plot %>% filter(measure=="EEG")
res = cor.test(data2use[,"param"], data2use[,"value"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use[, "param"] and data2use[, "value"]
## t = 5.1283, df = 9, p-value = 0.000621
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5459841 0.9639247
## sample estimates:
##       cor 
## 0.8631563
res$p.value
## [1] 0.0006209952
# ==============================================================================
# out of sample prediction
trn_eeg_data = hm4di_pow_data %>% filter(dset=="training") %>% filter(measure=="EEG")
val_eeg_data = hm4di_pow_data %>% filter(dset=="validation") %>% filter(measure=="EEG")
oos_res_eeg = out_of_sample_prediction(trn_data = trn_eeg_data, 
                                   val_data = val_eeg_data, 
                                   dv_var="fr_tot", 
                                   pred_var="value")
print(sprintf("EEG R-sq = %f",oos_res_eeg$rsq))
## [1] "EEG R-sq = 0.940921"
trn_lfp_data = hm4di_pow_data %>% filter(dset=="training") %>% filter(measure=="LFP")
val_lfp_data = hm4di_pow_data %>% filter(dset=="validation") %>% filter(measure=="LFP")
oos_res_lfp = out_of_sample_prediction(trn_data = trn_lfp_data, 
                                   val_data = val_lfp_data, 
                                   dv_var="fr_tot", 
                                   pred_var="value")
print(sprintf("LFP R-sq = %f",oos_res_lfp$rsq))
## [1] "LFP R-sq = 0.964205"
# make plot of actual vs predicted
data2plot = rbind(oos_res_eeg$val_data, oos_res_lfp$val_data)

vlineval = NULL
yLabel = "Predicted Firing Rate (Hz)"
xLabel = "Actual Firing Rate (Hz)"
cols2use = dtype_cols2use
plot_diag_line = TRUE

pow_hm4di_El_oos_plot = make_scatterplot(data2plot = data2plot, 
                            x_var = "fr_tot",
                            y_var = "pred",
                            strat_var = "measure", 
                            vlineval = vlineval, 
                            xlabel2use = xLabel, 
                            ylabel2use = yLabel, 
                            title2use = title2use, 
                            cols2use = cols2use, 
                            plot_diag_line = plot_diag_line)
pow_hm4di_El_oos_plot

Plots

fig1 = (FR_g_plot + FR_camk_El_plot + FR_hm4di_El_plot) / 
  (H_g_plot + H_camk_El_plot + H_hm4di_El_plot) / 
  (H_g_oos_plot + H_camk_El_oos_plot + H_hm4di_El_oos_plot)
ggsave(filename = file.path(plotdir, "Figure01.pdf"), plot = fig1, width = 12, height = 12)
fig1

supp_fig1a = (slope_g_plot + slope_camk_El_plot + slope_hm4di_El_plot) / 
  (slope_g_oos_plot + slope_camk_El_oos_plot + slope_hm4di_El_oos_plot)
ggsave(filename = file.path(plotdir, "SuppFig01a.pdf"), plot = supp_fig1a, width = 12, height = 8)
supp_fig1a

supp_fig1b = (pow_g_plot + pow_camk_El_plot + pow_hm4di_El_plot) / 
  (pow_g_oos_plot + pow_camk_El_oos_plot + pow_hm4di_El_oos_plot)
ggsave(filename = file.path(plotdir, "SuppFig01b.pdf"), plot = supp_fig1b, width = 12, height = 8)
supp_fig1b